[1.] Imágenes satelitales

1.1 Motivación

Night lights and economic growth

Asian Financial Crisis: Java, Indonesia

Tomado de: Measuring Economic Growth from Outer Space

Genocide Event: Rwanda

Tomado de: Measuring Economic Growth from Outer Space

Puede acceder al documento aquí: Measuring Economic Growth from Outer Space

Measuring the size and growth of cities using nighttime light

Puede acceder al documento aquí: Measuring the size and growth of cities using nighttime light

Housing Vacancy Rate (VHR) y Luces Nocturnas

Puede acceder al documento aquí: An estimation of housing vacancy rate using NPP-VIIRS night-time light data and OpenStreetMap data

CAUSES OF SPRAWL: A PORTRAIT FROM SPACE

Urban Land in Miami, FL

Tomado de: CAUSES OF SPRAWL: A PORTRAIT FROM SPACE

Puede acceder al documento aquí: CAUSES OF SPRAWL: A PORTRAIT FROM SPACE

1.2 ¿Datos raster?

¿Qué es un raster?

Resolución

Bandas de un raster

1.3 Fuentes:

1.3.1 Night lights

DMSPL (1992-2013)

Puede acceder a los datos anuales aquí 1992-2103

VIIRS (2012-Currently)

  • Puede acceder a la web del proyecto aquí: link

  • Puede acceder a los datos mensuales aquí: link

  • Puede acceder a los datos diarios aquí: link

Harmonized (1992-2020)

Puede acceder a los datos y al paper aquí: link

1.3.2 Land cover

  • Puede acceder a la web del proyecto aquí: link

  • Puede acceder a los datos aquí: link

  • Puede acceder a la documentación del proyecto aquí: link ## Dictionary

1.3.3 NASA

Puede acceder a datos de lluvias, temperatura, contaminación, vientos y otros en la página oficial del proyecto GES DISC de la NASA aquí: https://disc.gsfc.nasa.gov.

[2.] Manipular datos raster en R

0.0 Configuración inicial

Para replicar las secciones de esta clase, primero debe descargar el siguiente proyecto de R y abrir el archivo clase-05.Rproj.

Instalar/llamar las librerías de la clase

## Instalar/llamar las librerías de la clase
require(pacman) 
p_load(tidyverse,rio,skimr,viridis,osmdata,
       raster,stars, ## datos raster
       ggmap, ## get_stamenmap
       spatstat, ## analis de clusters
       sf, ## Leer/escribir/manipular datos espaciales
       leaflet) ## Visualizaciones dinámicas

## solucionar conflictos de funciones
select = dplyr::select

2.1 Leer un raster

methods(class = "stars")
##  [1] [                 [[<-              [<-               %in%             
##  [5] $<-               adrop             aggregate         aperm            
##  [9] as_tibble         as.data.frame     as.owin           c                
## [13] coerce            contour           cut               dim              
## [17] dimnames          dimnames<-        droplevels        filter           
## [21] hist              image             initialize        is.na            
## [25] Math              merge             mutate            Ops              
## [29] plot              predict           print             pull             
## [33] replace_na        select            show              slice            
## [37] slotsFromS3       split             st_apply          st_area          
## [41] st_as_sf          st_as_sfc         st_as_stars       st_bbox          
## [45] st_coordinates    st_crop           st_crs            st_crs<-         
## [49] st_dimensions     st_dimensions<-   st_extract        st_geometry      
## [53] st_interpolate_aw st_intersects     st_join           st_mosaic        
## [57] st_normalize      st_redimension    st_sample         st_set_bbox      
## [61] st_transform_proj st_transform      transmute         write_stars      
## see '?methods' for accessing help and source code
## importar raster de luces: raster
luces_r = raster('input/raster/night_light_202301.tif')
luces_r
## class      : RasterLayer 
## dimensions : 2990, 2919, 8727810  (nrow, ncol, ncell)
## resolution : 0.004166667, 0.004166667  (x, y)
## extent     : -79.01042, -66.84792, 0.002082733, 12.46042  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : night_light_202301.tif 
## names      : night_light_202301 
## values     : -0.64, 2208.49  (min, max)
## importar raster de luces: stars
luces_s = read_stars("input/raster/night_light_202301.tif")
luces_s
## stars object with 2 dimensions and 1 attribute
## attribute(s), summary of first 1e+05 cells:
##                         Min. 1st Qu. Median      Mean 3rd Qu. Max.  NA's
## night_light_202301.tif  0.07    0.16   0.18 0.1810952     0.2 0.53 97553
## dimension(s):
##   from   to   offset       delta refsys point values x/y
## x    1 2919 -79.0104  0.00416667 WGS 84 FALSE   NULL [x]
## y    1 2990  12.4604 -0.00416667 WGS 84 FALSE   NULL [y]

2.2 Atributos de un raster

## resolucion
0.00416667*110000 
## [1] 458.3337
## geometria
st_bbox(luces_s)
##          xmin          ymin          xmax          ymax 
## -79.010415858   0.002082733 -66.847915761  12.460416166
st_crs(luces_s)
## Coordinate Reference System:
##   User input: WGS 84 
##   wkt:
## GEOGCRS["WGS 84",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     ID["EPSG",4326]]
st_dimensions(luces_s)
##   from   to   offset       delta refsys point values x/y
## x    1 2919 -79.0104  0.00416667 WGS 84 FALSE   NULL [x]
## y    1 2990  12.4604 -0.00416667 WGS 84 FALSE   NULL [y]
## atributos
names(luces_s)
## [1] "night_light_202301.tif"
names(luces_s) = "date_202301"

## valores del raster
luces_s[[1]] %>% max(na.rm = T)
## [1] 2208.49
luces_s[[1]] %>% min(na.rm = T)
## [1] -0.64
luces_s[[1]] %>% as.vector() %>% summary() 
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##      -1       0       0       0       0    2208 4031249
luces_s[[1]][is.na(luces_s[[1]])==T] %>% head()# Reemplazar NA's
## [1] NA NA NA NA NA NA
luces_s[[1]][2000:2010,2000:2010]
##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
##  [1,] 0.15 0.12 0.09 0.17 0.14 0.12 0.18 0.20 0.24  0.28  0.14
##  [2,] 0.15 0.12 0.07 0.16 0.14 0.13 0.12 0.16 0.16  0.25  0.14
##  [3,] 0.06 0.10 0.13 0.20 0.15 0.11 0.10 0.09 0.10  0.17  0.12
##  [4,] 0.07 0.13 0.18 0.16 0.11 0.20 0.21 0.08 0.12  0.14  0.12
##  [5,] 0.08 0.13 0.16 0.16 0.15 0.22 0.21 0.19 0.15  0.13  0.11
##  [6,] 0.12 0.13 0.16 0.15 0.13 0.18 0.19 0.17 0.15  0.10  0.13
##  [7,] 0.12 0.15 0.19 0.24 0.19 0.22 0.15 0.11 0.14  0.13  0.13
##  [8,] 0.16 0.11 0.08 0.18 0.21 0.19 0.11 0.15 0.13  0.19  0.19
##  [9,] 0.13 0.12 0.19 0.16 0.13 0.15 0.19 0.15 0.09  0.13  0.17
## [10,] 0.15 0.15 0.16 0.16 0.11 0.06 0.07 0.14 0.11  0.11  0.18
## [11,] 0.11 0.13 0.16 0.17 0.14 0.08 0.06 0.17 0.18  0.17  0.21
luces_s[[1]][2000:2010,2000:2010] %>% table() # Sustraer una parte de la matriz
## .
## 0.0599999986588955 0.0700000002980232 0.0799999982118607 0.0900000035762787 
##                  3                  3                  4                  3 
##  0.100000001490116  0.109999999403954  0.119999997317791  0.129999995231628 
##                  4                 10                 10                 15 
##  0.140000000596046  0.150000005960464  0.159999996423721  0.170000001788139 
##                  8                 14                 12                  7 
##  0.180000007152557  0.189999997615814  0.200000002980232  0.209999993443489 
##                  6                  9                  3                  4 
##  0.219999998807907  0.239999994635582               0.25  0.280000001192093 
##                  2                  2                  1                  1
## puedo reproyectar un raster?
st_crs(luces_s)
## Coordinate Reference System:
##   User input: WGS 84 
##   wkt:
## GEOGCRS["WGS 84",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     ID["EPSG",4326]]
luces_new_crs = st_transform(luces_s,crs=4126)
luces_s[[1]][2000:2010,2000:2010] # no se alteran las geometrias
##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
##  [1,] 0.15 0.12 0.09 0.17 0.14 0.12 0.18 0.20 0.24  0.28  0.14
##  [2,] 0.15 0.12 0.07 0.16 0.14 0.13 0.12 0.16 0.16  0.25  0.14
##  [3,] 0.06 0.10 0.13 0.20 0.15 0.11 0.10 0.09 0.10  0.17  0.12
##  [4,] 0.07 0.13 0.18 0.16 0.11 0.20 0.21 0.08 0.12  0.14  0.12
##  [5,] 0.08 0.13 0.16 0.16 0.15 0.22 0.21 0.19 0.15  0.13  0.11
##  [6,] 0.12 0.13 0.16 0.15 0.13 0.18 0.19 0.17 0.15  0.10  0.13
##  [7,] 0.12 0.15 0.19 0.24 0.19 0.22 0.15 0.11 0.14  0.13  0.13
##  [8,] 0.16 0.11 0.08 0.18 0.21 0.19 0.11 0.15 0.13  0.19  0.19
##  [9,] 0.13 0.12 0.19 0.16 0.13 0.15 0.19 0.15 0.09  0.13  0.17
## [10,] 0.15 0.15 0.16 0.16 0.11 0.06 0.07 0.14 0.11  0.11  0.18
## [11,] 0.11 0.13 0.16 0.17 0.14 0.08 0.06 0.17 0.18  0.17  0.21
luces_new_crs[[1]][2000:2010,2000:2010] # no se alteran las geometrias
##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
##  [1,] 0.15 0.12 0.09 0.17 0.14 0.12 0.18 0.20 0.24  0.28  0.14
##  [2,] 0.15 0.12 0.07 0.16 0.14 0.13 0.12 0.16 0.16  0.25  0.14
##  [3,] 0.06 0.10 0.13 0.20 0.15 0.11 0.10 0.09 0.10  0.17  0.12
##  [4,] 0.07 0.13 0.18 0.16 0.11 0.20 0.21 0.08 0.12  0.14  0.12
##  [5,] 0.08 0.13 0.16 0.16 0.15 0.22 0.21 0.19 0.15  0.13  0.11
##  [6,] 0.12 0.13 0.16 0.15 0.13 0.18 0.19 0.17 0.15  0.10  0.13
##  [7,] 0.12 0.15 0.19 0.24 0.19 0.22 0.15 0.11 0.14  0.13  0.13
##  [8,] 0.16 0.11 0.08 0.18 0.21 0.19 0.11 0.15 0.13  0.19  0.19
##  [9,] 0.13 0.12 0.19 0.16 0.13 0.15 0.19 0.15 0.09  0.13  0.17
## [10,] 0.15 0.15 0.16 0.16 0.11 0.06 0.07 0.14 0.11  0.11  0.18
## [11,] 0.11 0.13 0.16 0.17 0.14 0.08 0.06 0.17 0.18  0.17  0.21
## plot data
plot(luces_s)
## downsample set to c(7,7)

2.3 Filtrar usando la gemoetría

## download boundary
quilla <- getbb(place_name = "Barranquilla, Colombia", 
                featuretype = "boundary:administrative", 
                format_out = "sf_polygon") %>% .[1,]
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
leaflet() %>%
addTiles() %>%
addPolygons(data=quilla)
## cliping
l_quilla_1 = st_crop(x = luces_s , y = quilla) # crop luces de Colombia con polygono de Medellin

l_quilla_1[[1]] %>% t()
##       [,1]  [,2]  [,3]  [,4]  [,5]   [,6]   [,7]   [,8]   [,9] [,10]  [,11]
##  [1,]   NA    NA    NA    NA    NA     NA     NA     NA     NA    NA     NA
##  [2,]   NA    NA    NA    NA    NA  32.05  72.30     NA     NA    NA     NA
##  [3,]   NA    NA    NA    NA    NA     NA  94.43  92.54  72.77    NA     NA
##  [4,]   NA    NA    NA    NA    NA  95.34  76.51  56.01  59.52 47.77     NA
##  [5,]   NA    NA    NA    NA    NA  62.60  51.57  28.10  32.52 63.07  76.49
##  [6,]   NA    NA    NA    NA 15.24  21.01  27.16  46.55  50.28 71.32  83.99
##  [7,]   NA    NA    NA    NA 40.83  47.05  64.18  76.44  77.28 78.88  65.12
##  [8,]   NA    NA    NA 75.98 71.70  85.57 105.79  99.95  84.11 79.16  66.23
##  [9,]   NA    NA    NA 82.99 93.33 101.92 111.54 102.47  76.40 81.48  79.32
## [10,]   NA    NA    NA 76.66 85.32  91.35  99.35  96.30  91.31 86.57  86.37
## [11,]   NA    NA 34.13 62.03 71.77  79.13  87.62  95.27 100.51 87.12  86.48
## [12,]   NA    NA    NA 35.22 45.27  62.79  76.90  89.58  96.36 92.16  90.84
## [13,]   NA    NA    NA 44.63 33.76  42.14  58.06  68.38  83.42 87.67  95.82
## [14,]   NA    NA    NA 43.40 22.66  35.22  49.98  60.20  68.46 79.86 100.32
## [15,]   NA    NA    NA    NA 42.55  45.05  52.81  58.82  61.30 65.44  69.87
## [16,]   NA    NA    NA    NA 62.58  64.82  55.21  55.67  55.18 55.81  61.82
## [17,]   NA 36.12 36.38 49.91 70.47  69.25  54.57  48.39  55.60 57.99  59.12
## [18,]   NA 39.68 49.92 55.94 69.37  62.20  45.75  43.38  48.67 58.43  58.08
## [19,]   NA    NA 51.42 59.80 58.90  48.16  43.53  42.76  48.42 56.79  54.78
## [20,]   NA    NA 47.31 63.17 52.35  46.63  45.12  50.09  56.28 57.06  53.92
## [21,]   NA    NA 66.26 65.72 53.72  47.85  49.67  57.15  60.50 53.60  50.88
## [22,]   NA    NA 70.88    NA 56.34  52.74  53.29  67.56  65.18 53.31  55.93
## [23,]   NA    NA    NA    NA 57.27  55.93  56.22  63.81  57.03 52.76  59.40
## [24,]   NA    NA    NA    NA 54.40  56.93  57.73  57.26  51.44 49.65  50.65
## [25,]   NA    NA    NA    NA 60.24  69.57  64.83  56.72  47.53 44.60  44.62
## [26,]   NA    NA    NA    NA 54.39  69.89  67.61  56.60  48.08 45.28  47.42
## [27,]   NA    NA    NA    NA 47.97  61.82  63.79  62.72  59.50 49.69  52.51
## [28,]   NA    NA    NA    NA 25.88  43.63  49.79  63.71  60.43 52.18  59.41
## [29,]   NA    NA    NA    NA    NA  47.43  54.56  61.69  58.86 60.99  69.98
## [30,]   NA    NA    NA    NA    NA     NA  58.19  64.35  67.04 71.71  82.92
## [31,]   NA    NA    NA    NA    NA     NA     NA     NA     NA    NA     NA
##       [,12]  [,13] [,14] [,15] [,16] [,17]  [,18] [,19] [,20] [,21] [,22] [,23]
##  [1,]    NA     NA    NA    NA    NA    NA     NA    NA    NA    NA    NA    NA
##  [2,]    NA     NA    NA    NA    NA    NA     NA    NA    NA    NA    NA    NA
##  [3,]    NA     NA    NA    NA    NA    NA     NA    NA    NA    NA    NA    NA
##  [4,]    NA     NA    NA    NA    NA    NA     NA    NA    NA    NA    NA    NA
##  [5,] 85.56     NA    NA    NA    NA    NA     NA    NA    NA    NA    NA    NA
##  [6,] 85.30  63.33    NA    NA    NA    NA     NA    NA    NA    NA    NA    NA
##  [7,] 60.87  54.69 53.39    NA    NA    NA     NA    NA    NA    NA    NA    NA
##  [8,] 53.39  47.17 49.70 51.41    NA    NA     NA    NA    NA    NA    NA    NA
##  [9,] 66.51  56.38 59.51 62.07 46.97    NA     NA    NA    NA    NA    NA    NA
## [10,] 84.80  69.04 69.33 66.04 51.06 27.09     NA    NA    NA    NA    NA    NA
## [11,] 85.51  77.50 65.75 62.99 51.63 51.20  36.54    NA    NA    NA    NA    NA
## [12,] 87.09  82.91 75.94 66.18 67.95 74.07  70.44 46.73    NA    NA    NA    NA
## [13,] 90.66  76.26 78.83 78.96 75.55 81.83  70.24 36.41 18.83    NA    NA    NA
## [14,] 91.54  73.72 75.44 93.24 92.60 77.20  65.92 41.21 23.91 22.20    NA    NA
## [15,] 71.60  71.90 75.94 88.63 98.86 86.89  71.64 60.35 49.21 42.73    NA    NA
## [16,] 66.07  71.64 73.51 78.55 92.67 85.76  82.93 82.86 69.00 61.91 51.15    NA
## [17,] 65.07  69.80 67.73 73.24 83.59 81.23  83.38 95.06 82.41 64.15 52.33    NA
## [18,] 60.75  62.08 58.91 66.09 72.25 70.98  76.85 87.11 81.52 44.85 37.96    NA
## [19,] 54.65  60.34 60.78 67.11 64.15 64.04  64.97 68.56 65.19 51.40 72.00    NA
## [20,] 57.21  72.60 79.05 76.00 68.02 64.23  63.44 70.97 66.07 56.48    NA    NA
## [21,] 54.49  66.49 70.52 66.77 71.46 68.73  68.48 75.41 76.00 55.55    NA    NA
## [22,] 60.69  63.25 60.49 60.49 71.85 89.37  73.30 72.00 73.48 63.38 50.81    NA
## [23,] 61.73  75.28 80.32 70.60 82.59 84.18  70.54 71.74 76.27 80.48 83.79    NA
## [24,] 51.92  69.87 86.02 78.78 92.53 87.69  74.86 83.80 94.17 93.37 85.02    NA
## [25,] 49.16  59.19 73.68 80.36 92.06 99.48  96.92 87.09 91.57 88.35 76.77    NA
## [26,] 49.30  67.65 71.78 69.61 84.80 96.47 111.94    NA    NA    NA    NA    NA
## [27,] 58.76  80.24 76.63 67.45 72.62 79.86  97.51    NA    NA    NA    NA    NA
## [28,] 72.81  93.18 94.09 87.39    NA    NA     NA    NA    NA    NA    NA    NA
## [29,] 79.17  96.70 97.20    NA    NA    NA     NA    NA    NA    NA    NA    NA
## [30,] 91.30 119.18    NA    NA    NA    NA     NA    NA    NA    NA    NA    NA
## [31,]    NA     NA    NA    NA    NA    NA     NA    NA    NA    NA    NA    NA
## Plot data
ggplot() + geom_stars(data=l_quilla_1 , aes(y=y,x=x,fill=date_202301)) + # plot raster
scale_fill_viridis(option="A" , na.value='white') +
geom_sf(data=quilla , fill=NA , col="green") + theme_bw() 

[3.] Aplicación: Actividad económica en Barranquilla

3.1 Importar datos 201301

## load data
l_quilla_0 = read_stars("input/raster/night_light_201301.tif") %>% st_crop(quilla)
names(l_quilla_0) = "date_201301"

ggplot() + geom_stars(data=l_quilla_0 , aes(y=y,x=x,fill=date_201301)) + # plot raster
scale_fill_viridis(option="A" , na.value='white') +
geom_sf(data=quilla , fill=NA , col="green") + theme_bw() 

## stack rasters: adicionar bandas
l_quilla = c(l_quilla_0,l_quilla_1)
l_quilla
## stars object with 2 dimensions and 2 attributes
## attribute(s):
##               Min. 1st Qu. Median     Mean 3rd Qu.   Max. NA's
## date_201301   8.70   55.10  65.90 66.55410   78.99 135.63  308
## date_202301  15.24   54.49  65.72 66.38741   78.83 119.18  308
## dimension(s):
##   from   to   offset       delta refsys point values x/y
## x  999 1021 -79.0104  0.00416667 WGS 84 FALSE   NULL [x]
## y  340  370  12.4604 -0.00416667 WGS 84 FALSE   NULL [y]

3.2 Raster a datos vectoriales

## st_as_sf
puntos_quilla = st_as_sf(x = l_quilla, as_points = T, na.rm = T) # raster to sf (points)
puntos_quilla
## Simple feature collection with 405 features and 2 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -74.85 ymin: 10.92083 xmax: -74.75833 ymax: 11.04583
## Geodetic CRS:  WGS 84
## First 10 features:
##    date_201301 date_202301                   geometry
## 1        34.77       32.05 POINT (-74.82917 11.04167)
## 2        61.45       72.30   POINT (-74.825 11.04167)
## 3        73.49       94.43    POINT (-74.825 11.0375)
## 4        78.52       92.54  POINT (-74.82083 11.0375)
## 5        57.54       72.77  POINT (-74.81667 11.0375)
## 6        32.77       95.34 POINT (-74.82917 11.03333)
## 7        28.12       76.51   POINT (-74.825 11.03333)
## 8        27.86       56.01 POINT (-74.82083 11.03333)
## 9        36.19       59.52 POINT (-74.81667 11.03333)
## 10       45.98       47.77  POINT (-74.8125 11.03333)
poly_quilla = st_as_sf(x = l_quilla, as_points = F, na.rm = T) # raster to sf (polygons)
poly_quilla
## Simple feature collection with 405 features and 2 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -74.85208 ymin: 10.91875 xmax: -74.75625 ymax: 11.04792
## Geodetic CRS:  WGS 84
## First 10 features:
##    date_201301 date_202301                       geometry
## 1        34.77       32.05 POLYGON ((-74.83125 11.0437...
## 2        61.45       72.30 POLYGON ((-74.82708 11.0437...
## 3        73.49       94.43 POLYGON ((-74.82708 11.0395...
## 4        78.52       92.54 POLYGON ((-74.82292 11.0395...
## 5        57.54       72.77 POLYGON ((-74.81875 11.0395...
## 6        32.77       95.34 POLYGON ((-74.83125 11.0354...
## 7        28.12       76.51 POLYGON ((-74.82708 11.0354...
## 8        27.86       56.01 POLYGON ((-74.82292 11.0354...
## 9        36.19       59.52 POLYGON ((-74.81875 11.0354...
## 10       45.98       47.77 POLYGON ((-74.81458 11.0354...
## Plot data
ggplot() + 
geom_sf(data = puntos_quilla , aes(col=date_202301))  + 
scale_color_viridis(option="A" , na.value='white') +
geom_sf(data=quilla , fill=NA , col="green") + theme_test() 

ggplot() + 
geom_sf(data = poly_quilla , aes(fill=date_202301),col=NA)  + 
scale_fill_viridis(option="A" , na.value='white') +
geom_sf(data=quilla , fill=NA , col="green") + theme_test()

3.3 Asignar el valor del pixel a un lugar

## get polygon
puerto <- getbb(place_name = "Puerto de Barranquilla - Sociedad Portuaria, Colombia", 
                featuretype = "barrier:wall", 
                format_out = "sf_polygon")
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
leaflet() %>%
addTiles() %>%
addPolygons(data=puerto)
## asignar luminosidad
l_puerto = st_join(puerto,poly_quilla)

## variacion promedio 
df = l_puerto
st_geometry(df) = NULL

df %>% 
summarise(pre=mean(date_201301,na.rm=T) , 
          post=mean(date_202301,na.rm=T)) %>%
mutate(ratio=post/pre-1)
##    pre  post      ratio
## 1 61.6 54.46 -0.1159091

[4.] Aplicación: Cambio de cobertura de suelo?

## read raster
land_cover = c(read_stars("input/raster/discreta_2015.tif"),
               read_stars("input/raster/discreta_2019.tif"))
0.000992063*111000 # resolucion
## [1] 110.119
## read polygon
quilla_buffer <- st_buffer(x = quilla,dist = 2000)

leaflet() %>%
addTiles() %>%
addPolygons(data=quilla_buffer)
quilla_border <- st_difference(x = quilla_buffer , y = quilla)

leaflet() %>%
addTiles() %>%
addPolygons(data=quilla_border)
## cliping raster
quilla_cover = land_cover %>% st_crop(quilla_buffer)
plot(quilla_cover) ## view data

## cliping raster
border = land_cover %>% st_crop(quilla_border)
plot(border) ## view data

## raster to sf
border_sf = st_as_sf(x=border, as_points = F, na.rm = T) 

## cuantos pixeles se urbanizaron?
border_sf = border_sf %>% 
            rename(c_2015=discreta_2015.tif,c_2019=discreta_2019.tif)
border_sf = border_sf %>% 
            mutate(build = ifelse(c_2015!=50 & c_2019==50,1,0))
table(border_sf$build)
## 
##    0    1 
## 7970   52

[5.] Aplicación: Clusteres espaciales

## restaurantes
points <- opq(bbox = getbb("Bogota Colombia")) %>%
          add_osm_feature(key = "amenity", value = "restaurant") %>%
          osmdata_sf() %>% .$osm_points %>% select(osm_id,name) %>% mutate(restaurant=1)

leaflet() %>% addTiles() %>% addCircles(data=points)
## download boundary
city <- getbb(place_name = "Bogota, Colombia", 
        featuretype = "boundary:administrative", 
        format_out = "sf_polygon") %>% .[[2]]
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
leaflet() %>% addTiles() %>% addPolygons(data=city)
## afine transformations
points = st_transform(points,3857)
city = st_transform(city,3857)

## sf to ppp
points_p <- as.ppp(X = points["restaurant"])
Window(points_p) <- as.owin(city)
plot(points_p)

#--------------#
## summary
summary(points_p) # las unidades son en metros
## Marked planar point pattern:  6751 points
## Average intensity 1.697305e-05 points per square unit
## 
## Coordinates are given to 1 decimal place
## i.e. rounded to the nearest multiple of 0.1 units
## 
## marks are numeric, of type 'double'
## Summary:
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##       1       1       1       1       1       1 
## 
## Window: polygonal boundary
## 5 separate polygons (no holes)
##            vertices        area relative.area
## polygon 1      7588 395800000.0      0.995000
## polygon 2       276   1524690.0      0.003830
## polygon 3        72    235244.0      0.000591
## polygon 4        45     75844.4      0.000191
## polygon 5        81    112582.0      0.000283
## enclosing rectangle: [-8262524, -8238784] x [498234.9, 538665.2] units
##                      (23740 x 40430 units)
## Window area = 397748000 square units
## Fraction of frame area: 0.414
intensity(points_p)
## [1] 1.697305e-05
## reescalar
points_p = rescale(X = points_p , s = 100 , unitname = "100-metros")
unitname(points_p)
## 100-metros / 100-metros
summary(points_p)
## Marked planar point pattern:  6751 points
## Average intensity 0.1697305 points per square 100-metros
## 
## Coordinates are given to 3 decimal places
## i.e. rounded to the nearest multiple of 0.001 100-metros
## 
## marks are numeric, of type 'double'
## Summary:
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##       1       1       1       1       1       1 
## 
## Window: polygonal boundary
## 5 separate polygons (no holes)
##            vertices        area relative.area
## polygon 1      7588 39580.00000      0.995000
## polygon 2       276   152.46900      0.003830
## polygon 3        72    23.52440      0.000591
## polygon 4        45     7.58444      0.000191
## polygon 5        81    11.25820      0.000283
## enclosing rectangle: [-82625.24, -82387.84] x [4982.349, 5386.652] 100-metros
##                      (237.4 x 404.3 100-metros)
## Window area = 39774.8 square 100-metros
## Unit of length: 1 100-metros
## Fraction of frame area: 0.414
intensity(points_p)
## [1] 0.1697305
## estimated standard error for intensity
sqrt(intensity(points_p)/summary(points_p)$window$area)
## [1] 0.002065741
## plot 
plot(density(points_p,5)) ## plot density

contour(density(points_p,bw.ppl(points_p)), axes = FALSE) ## contour density

#--------------#
## quadrants
points_q = quadratcount(points_p, nx = 4, ny = 6)
points_q
## tile
## Tile row 1, col 3 Tile row 1, col 4 Tile row 2, col 2 Tile row 2, col 3 
##                 0                20               110               335 
## Tile row 2, col 4 Tile row 3, col 1 Tile row 3, col 2 Tile row 3, col 3 
##               759                 0               362               835 
## Tile row 3, col 4 Tile row 4, col 1 Tile row 4, col 2 Tile row 4, col 3 
##               806                54               650              2031 
## Tile row 4, col 4 Tile row 5, col 1 Tile row 5, col 2 Tile row 5, col 3 
##               237                13               187               341 
## Tile row 5, col 4 Tile row 6, col 2 Tile row 6, col 3 
##                 0                 0                11
plot(points_q, cex = 2 , col="red")

Referencias

  • Lovelace, R., Nowosad, J., & Muenchow, J. (2019). Geocomputation with R. [Ver aquí]

    • Cap. 4: Spatial Data Operations
    • Cap. 5: Geometry Operations
    • Cap. 6: Reprojecting geographic data
    • Cap. 11: Statistical learning
  • Bivand, R. S., Pebesma, E. J., Gómez-Rubio, V., & Pebesma, E. J. (2013). Applied spatial data analysis with R. [Ver aquí]

    • Cap. 4: Spatial Data Import and Export
    • Cap. 7: Spatial Point Pattern Analysis
    • Cap. 8: Interpolation and Geostatistics